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Abstract 

An efficient geometric integrator is proposed for solving the perturbed Kepler mo- 
tion. This method is stable and accurate over long integration time, which makes it 
appropriate for treating problems in astrophysics, like solar system simulations, and 
atomic and molecular physics, like classical simulations of highly excited atoms in 
external fields. The key idea is to decompose the hamiltonian in solvable parts and 
propagate the system according to each term. Two case studies, the Kepler atom 
in an uniform field and in a monochromatic field, are presented and the errors are 
analyzed. 
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1 Introduction 



The motion of a body under a force which depends inverse quadratically with 
the distance is known as the Kepler problem and explains both the planetary 
dynamics and the electronic structure of atoms. The problem is completely 
solvable and is a textbook example of constructing action-angle variables for 
a nontrivial system. The stable computation of trajectories of perturbed Ke- 
pler systems is important for astronomical applications, such as solar system 
studies, and in the classical and semi-classical studies of atomic systems, par- 
ticularly atoms in highly excited (Rydberg) states. Once the initial conditions 
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are assigned, the resulting system of ordinary differential equations can be nu- 
merically solved using standard, good-for-all, algorithms such as Runge-Kutta 
or Gear integrators. Although such methods have a good control over the local 
error, their global error usually grows exponentially. This is why the traditional 
numerical methods are unsuitable in treating physical problems which require 
long time integration. The simulation of the solar system for more than 10 9 
years [1] is an example where special integration methods are needed [2]. Due 
to the weak electromagnetic coupling (of the order of fine structure constant 
cubed ps 1/137 3 ) the simulation of the interaction of a Rydberg atom with 
radiation also requires following the electron trajectory for a long time. 

Geometric integration is a relatively new branch of the Computational Math- 
ematics. It studies algorithms and discretization methods that respect the 
underlying geometry and qualitative structure of the the problems it aims to 
solve. The principle is that if more specific information about the problem is 
explicitly included into the solver, then the solution will be more accurate and 
stable than those produced by generic methods. As a simple illustration, con- 
sider the harmonic oscillator problem du/dt = v, dv/dt = —u. This system has 
the property that (a) its solutions are all periodic and, (b) bounded and, (c) 
that u 2 +v 2 is a conserved quantity of the evolution. A forward Euler discretiza- 
tion of this system with time step h gives u n+ i = u n + hv n , v n+ i = v n — hu n . It 
is easy to see that u^ +1 +v 2 +1 = (1 + h 2 )(u^ + 0- While the correct solution 
is obtained in the limit of h — > 0, the numerical method, over any finite time 
interval, has lost all three qualitative features listed above! 

For a Hamiltonian system, a numerical method is symplectic if it produces an 
approximate solution and in the same time it explicitly preserves the underly- 
ing symplectic structure. Many symplectic methods also obey other qualities 
of the system such as symmetry and time reversal. Although the Hamiltonian 
itself is not conserved in general (for autonomous systems), it was demon- 
strated that the error grows only linearly in time. Comprehensive reviews of 
symplectic methods are presented in references [3] and [4]. 

Suppose that the Hamiltonian is separable and each term in the separation 
defines an integrable problem, then a symplectic solution is obtained by propa- 
gating the system with each term separately. The symplectic quality is evident 
since the composition of two symplectic maps is also symplectic and the solu- 
tion is exact in the limit of vanishing time step. The traditional way is to split 
the Hamiltonian as the kinetic energy T plus potential energy V. Evolution 
under the kinetic energy part alone describes the free motion of the system. 
The potential energy term usually does not depend on momentum coordinates 
and therefore the position coordinates remain unchanged during propagation 
with it. The system gets a "kick", a sudden change in momentum. It was 
proved (see [3], for example) that the leap-frog, or Stormer-Verlet, methods, 
very popular in Molecular Dynamics simulations, can be derived from &T + V 
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splitting and are therefore symplectic. Essential for a more general way of 
splitting the Hamiltonian is that the parts are not only integrable but also 
efficiently computable. 

This paper proposes a symplectic integration method for perturbed Kepler 
problems by splitting the Hamiltonian into a pure, exactly solvable, Kepler 
part and a perturbation part. Even when the perturbation is not small, the 
system can be propagated with accuracy using relatively large time steps. 
This method can be applied for all problems where such a decomposition 
makes sense. For a collision problem, such as electron scattering on an atom, 
this splitting is not practical when the projectile and target electrons come 
close to each other and interact strongly. In such a case the time step has to 
be chosen very small so that trajectories are straight lines (during one time 
step) and this method is no more efficient than the simpler T + V splitting. 
Although the algorithm is still correct, the advantage of being able to take 
large time steps is lost. 

Two case studies are presented to demonstrate the efficiency of the proposed 
method. Interaction of a Rydberg atom with an uniform external field is stud- 
ied in both cases. One case treats the interaction with a constant field (Stark 
effect) and the other with that of a monochromatic oscillatory field (laser 
interaction). 



2 Symplectic propagator 

The motion of a planet or of a Rydberg electron can be described by the 
classical equations of motion derived from the Hamiltonian 

H = H + V 

where H = p 2 /2 — 1/r (in an appropriate set of units) is the Kepler Hamilto- 
nian and V is the perturbation potential. The evolution in time of any function 
z defined over the phase space is given by the general equation 

dz 

- = {z,H} = D H (z) 

where { , } are the Poisson brackets and D H is an operator defined by the 
Hamiltonian H. Providing the initial condition z at to, the function z is 
obtained at time t as a solution of the above equation, given as z(t) = &At( z o)- 
Formally, the mapping $ can written as <3>a< — e AtE>H , in a form reminiscent 
of the quantum mechanical evolution operator. Of course, the exponential of 
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operator Dh has an exact meaning only if the solution of the Hamiltonian H is 
explicitly known. Apparently no real advantage is gained. However, solutions 
for H and V are separately known, and, for any At, the mapping associated 
with these solutions are e At£>H o and e AtDv . If the phase space operators D Ho 
and Dy commute (or their Poisson bracket {H , V} = 0) then the product of 
exponentials is equal to the exponential of their sum, and the problem has an 
exact solution: $At = e AtDH oe AtDv . This is not the case in general, and the 
following useful expansion can be derived from the Baker-Campbell-Hausdorff 
formula 



log(e AA / W A/2 ) = X(A + B)- — [AAB] - —[BAB] (1) 
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where the bracket notation refers to commutators: [XY] = [X, Y] = AB — BA, 
[XYZ] = [X, [Y, Z]], and so on. Operators A and B can be either D Ho or D v . 
The Poisson bracket also has the property of being a Lie brackets which means 
that [D X ,D Y ] = D {Y ^ } . 

The meaning of the expansion formula (1) is that the propagator obtained by 
successive application of e AA / 2 , e XB and t XA l 2 is a propagator of an equivalent 
Hamiltonian H equal to the original H = A + B, only in the limit of A — > 0. 
For a small enough time step A, the exact solution under H is expected to 
converge to a solution of H . This propagator is symplectic and also preserves 
the symmetry, time reversibility and first order invariants of the system. When 
the perturbation V does not depend on time then both Hamiltonians H and 
H are conserved. Therefore the global deviation from energy conservation is 
directly accessible from H - H. This is a valuable asset of the symplectic 
methods, since most good-for-all integration methods cannot easily predict 
the global behavior of their solutions, in general. 

Explicitly, from Eq. (1), on choosing A = D v and B = D Ho , one gets for the 
equivalent Hamiltonian 



H = H- ^{{H , V}, V + 2H } + 0(X 4 ). 



4 



Hence this propagator is second order in the time step A. The global error 
correct up the second order in A is then 

r( 2 ) = ^-iJ=-^[F 2 + 2F-F c + 2p-|-(p.F)) (2) 



where F = — <9V/<9r is the perturbation force and F c = — r/r 3 is the Coulomb 
force. For small values of r, the error is therefore dominated by A 2 F ■ r/r 3 . 
In contrast, the global error for a T + V splitting, obtained from Eq. (1) and 
using A = Dy-i/r and B = Dt, is 

r w = r»_^(p.. + 1! p.| ( p.pj) 



The truncation error for the T + V splitting is still second order in the time 
step, but is clearly inferior at small distance r, where it behaves as 1/r 4 

Formula (1) suggests a couple of ways to improve the performance of this 
method. First, if operator B is replaced by B' = B + \ 2 /2A [A+2B, [A, B]] and 
providing the the evolution under B' is explicitly known, then the A 3 term in 
expansion is removed and the error is now of the order of A 5 . The second way of 
accelerating the convergence is obtained by adapting the times step. Superior 
order methods are obtained by composing steps with appropriate variable time 
steps. Examples of such schemes are presented in the next section. 

A solution of the pertubed Kepler problem is therefore obtained by evolving 
the phase space point (r, p) of the system successively, under two elementary 
operations. The "drift stage" is the evolution under the Kepler Hamiltonian 
H , while the evolution under the perturbation Hamiltonian V is called the 
"kick stage" . 



2. 1 Drift stage 



Although a textbook example, an explicit solution of the Kepler problem is not 
entirely trivial. To state the problem, given the position r , and momentum 
Po at time to, one need to find positon r and momentum p at some other time 
t + At. This mapping denoted by D(At), gives the position and momentum 
along a Keplerian orbit, from an initial position and momentum, after a time 
At. 

The trajectory during the drift stage is a segment of an ellipse, a parabola or 
a hyperbola, depending whether the "local" energy w = H = p 2 /2 — 1/r is 
positive, zero or negative, respectively. The geometric size of the ellipse and 
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the orbital period also depend on the local energy. The angular momentum 
L = r x p and the Runge-Lenz A = (p 2 — 1/r) r — (rp) p vectors, are used to 
identify the orientation of the orbital plane in space. Although w, L and A can 
be calculated from the initial state vector (r ,p ), it is advantageous to keep 
these quantities in the state vector, alongside with position and momentum. 
This helps save a number of floating point operations. Besides it can also 
help avoid the accumulation of round-off errors that might creep up in the 
calculation. For instance the energy is sometimes obtained as a difference of 
two large numbers: the kinetic and potential energies. 

Having as input the thirteen dimensional state vector (r, p, w, L, A), the drift 
stage proceed as follows: 

(1) obtain the characteristic parameters of the orbit: semimajor axis a = 
l/2w, eccentricity e = ||A|| and orbital angular frequency uj = (2|u>|) 3//2 , 

(2) calculate the direction of the pericenter e.\ = A/e and the direction in 
the orbital plane perpendicular to it, e 2 = L x A/eL, 

(3) find the eccentric anomaly corresponding to the initial position uq = 
arctan(l — 2|iu|r, \j2\w\rp), 

(4) find the eccentric anomaly u after time At as u = Kepler(e, uo — e sin u + 
a; At) , where Kepler (e, M) is the solution of Kepler's equation u — e sin u — 
M = function of parameters e and M, 

(5) calculate the new position on the orbit corresponding to the new eccentric 
anomaly u as 

r = a(cosu — e) i\ + ay/l — e 2 sin-ue 2 

1 sinw /l — e 2 cosw 

P = ^\ ei + \\ e 2 , 

y/ a 1 — e cos u V a 1 — e cos u 

(6) energy, angular momentum and Runge-Lenz vectors are not modified 
during the drift stage. 

The steps above apply specifically to the case of negative energy (elliptic orbit). 
It is not difficult to generalize this procedure for the parabolic and hyperbolic 
motions. 

Up to round-off errors, the drift stage integrates exactly the orbit for any time 
step At, except for the solution of the transcendental Kepler's equation which 
has to be obtained approximately. However, this equation has been long and 
carefully studied, as Goldstein remarks [5]: 

Indeed, it can be claimed that the practical need to solve Kepler's equation 
to accuracies of a second of arc over the whole range of eccentricity fathered 
many of the developments in numerical mathematics in the eighteenth and 
nineteenth centuries. A few of the more than 100 methods of solution devel- 
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oped in the pre-computer era are considered in the exercises to this chapter. 
2.2 Kick stage 

In the kick stage the system evolves solely under the perturbation Hamiltonian 
V. This mapping is denoted by K(At). If V does not depend on momentum, 
then the position vector is a cyclic coordinate and does not change during this 
stage. The change in momentum can then be explicitly obtained as 



Where F is the perturbation force derived from the potential V. Quantities 
w, L and A in the state vector are updated from r, p and Ap, instead of 
calculating them directly from r and p'. This precludes the accumulation of 
round-off errors and increases the efficiency of the procedure. For example, 
energy is updated during the kick stage as 

w' = w + pAp + ^Ap 2 



2.3 Kepler solver 

The solution of the transcendental Kepler's equation is the most time con- 
suming part in the propagator. The traditional numerical scheme to obtain 
accurate solutions is to "guess" a good starting approximation and then refine 
it by using Newton-Raphson iterations until the desired accuracy is obtained. 
Each iteration involves evaluating trigonometric functions several times. Since 
each trigonometric function evaluation has a cost of at least several hundreds 
of microprocessor clocks, regardless if it is done "on-the-chip" or by a library 
call, the cost of solving the Kepler's equation can mount easily to thousands 
of clocks. It is clear that a long time integration, on the order of 10 8 time steps 
would require a more refined procedure. A table-driven procedure is proposed 
here, which trades memory space in favor of time. No trigonometric functions 
are calculated during iterations. 

The Kepler solver takes two branches, depending whether the orbit is elliptic 
or hyperbolic. For negative energy, or when < e < 1, the equation to solve 
is 

M-esinw-M = 0. (3) 



to+At 




Fdt 
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In the hyperbolic case, for positive energy, or for 1 < e < oo, Kepler's equation 
has the form 



Equation (3) is solved as follows: 

(1) using the fact that Eq. (3) is invariant to (M — > M + nn, u — > u + nit, 
e — > (— l) n e) and (M — > 7r — M, w — > n—u) transformations, the argument 
M of the equation can be mapped to the [0, ir) interval. The equation 
needs to be solved only for < e < 1 and < M < it. The solution for 
arbitrary M is obtained by adding nn. 

(2) for large eccentricity (e — > 1) and low M the solution of Eq. (3) has an 
essential singularity. This can be seen by seeking a solution as power se- 
ries in M. Equating the coefficients, order by order, one gets u(e, M) = 
M/(l - e) - M 4 e/6(1 - e) 4 + C(M 5 ). This series does not converge uni- 
formly when e — > 1. On the other hand when e = 1 and the Taylor 
expansion of sin is used one gets u(e = 1,M) = (6M) 1 / 3 + 0(M). A 
direct Newton-Raphson approach is started for the extreme cases when 
e-lsiMsiO, using (6M) 1 / 3 as initial guess. 

(3) for the "regular" cases, a grid is set up for the (0, ir) interval with points 
Uk = kir/n g , k = 0, . . . ,n g and a table of sine and cosine at the grid 
points is stored at the start of the program. From an initial grid index 
guess u = M, a Newton-Raphson iteration for indices is started using 
the following mapping: k — > k' = k + S[(M + es k — u k )/{\ — ec^)], where 
S is a function returning the integer truncation of a real number, and 
Sk and Cfc are the sine and cosine values at the grid point k. At the end 
of this process the solution is localized between two grid points with an 
accuracy no greater than the grid spacing n/n g . In this way, the process 
has both the quadratic convergence of Newton-Raphson's method and 
the stability of bisection method. 

(4) A last fifth-order Newton-Raphson refinement is obtained using the grid 
point closest to the solution obtained at the previous step. However, in- 
stead of searching for a u value which satisfies Eq. (3), it is more efficient 
to look for a solution in the unknown tan(«/2). This procedure deliv- 
ers then directly the pair (sinu, cosw) corresponding to the solution, and 
hence no trigonometric function need to be explicitly calculated during 
the procedure call! Assuming that the values of the function /o and of 
the first four derivatives /i, /2, fz and f\ are known at the grid point, 
then the following corrections are obtained 



e sinh u — u — M = 0. 



(4) 




fi + f h 



fo 



8 



fo 

fl + f/2 + f/3 

fo 

fl + f/2 + f/3 + 3/4 



to give the approximate solution tan-u/2 pa \J(1 — c k )/{\ + c k ) + <5 5 

The accuracy of this procedure is expected to be of the order (it /n g ) 5 . Indeed, 
in tests using n g = 1024 the error in finding both sinw and cosw where 10~ 14 — 
— 10~ 15 except for the range of e — > 1 and M — > where the error is no greater 
than 10" 12 . 



A similar approach is taken to solve Eq. (4). Although Kepler's equation for 
positive energy is not manifestly periodic, a mapping to the standard interval 
is obtained by observing that Eq. (4) can be written as eie M — e 2 e~ M — u— M = 
and that the following set of transformations 



u — > u + nD M — > M + e x 2 — > €i 2 e 



±nD 



leave this form invariant. In this way, the equation needs to be solved only in 
the standard interval [0, D) where D is arbitrary. For convenience, D is chosen 
D = 2.0. This interval is gridded and the exponential = exp-u^ is calculated 
at each grid point k. The Newton-Raphson grid iteration 



k —>k' = k + S[(M + u k - eie k - e 2 e k )/ (eie fe + e 2 e k - 1)] 



ends by identifying the grid point closest to the solution. The solution for e u is 
obtained from this grid point after a fifth order Newton-Raphson refinement. 
This subroutine returns an approximation for the (sinh-u, coshw). Precision 
levels similar to the elliptic case are obtained. 

When tables of sine, cosine and exponential are build at the set up of the pro- 
gram, the drift stage can then be performed without evaluation of any tran- 
scendental function, which brings significant improvements to the efficiency of 
the integrator. 
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3 Applications 



3.1 Kepler atom in uniform electric field 



In the absence of a perturbation, the kick stage reduces to identity and the 
dynamics is described only by the drift. The propagator was tested and the 
solution was seen to be practically exact (within the machine precision) even 
for extremely long time integration. A proper test of this symplectic integrator 
can only be done in the presence of a perturbation. The simplest perturbation, 
which is also completely integrable, is the constant and uniform force field. 

When a constant force F acts on the system, the Kepler orbit orbit starts to 
precess and change its eccentricity. The corresponding potential is V = — rF, 
such that the energy E = p 2 /2 — 1/r— rF is conserved. For forces F > (-E/2) 2 
the system can break away and ionize. 

The angular momentum and the Runge-Lenz vectors evolve in time according 
to 



dL „ d 

— = r x F — 

dt dt 



A--rx (rxFl 

2 v ; 



-FxL 

2 



If the orbit does not change appreciably over one period, then the average 
angular momentum and Runge-Lenz vectors obey equations 



because of the Pauli's replacement rule r — > — (3/2)A, and the fact that 
r x (r x F)|q 0. Within these assumptions, the slow changes in L and A 
are obtained, by solving the above system of equations, as 

(L) = cos(^Ft)L(0) + [1 - cos(^Ft)][F • L(0)]F + sin(^Ft)[F x A(0)](5) 

(A) = cos(-Ft)A(0) + [1 - cos(-Ft)][F • A(0)]F + sin(-Fi)[F x L(0)](6) 
2 2 2 



Both L and A therefore rotate about each other with a period of 4ir/3F. 
The simplest second order symplectic integrator (step 2) 
S(At) = K(At/2)D(At)K(At/2) 
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Fig. 1. The relative energy error is compared for long time integration using the 
step2 integrator and implicit Runge-Kutta of order 4 from GSL. The initial orbit 
has eccentricity e = 0.9 and energy w = —0.5. The electric field is 5.5 x 10 -3 oriented 
perpendicularly on the orbit plane. The fixed time step for step2 corresponds to 200 
steps per orbit. 

requires only one "drift" stage D and two "kick" operations K. The "drift" 
stage involves solving Kepler's equation and has a higher computational cost 
than the "kick" stage. An equivalent integrator DKD is not as efficient, be- 
cause it uses two "drift" stages. 

Figure 1 compares the performance of step2 integrator with a standard implicit 
Runge-Kutta method of order 4, with adaptive time step (rk4imp), from the 
free GNU Scientific Library (GSL) [6]. The initial orbit has eccentricity 0.9, 
energy -0.5 and period 2ir in the chosen units. An uniform electric field of 
magnitude 5.5 x 10~ 3 is applied along a direction perpendicular to the orbital 
plane. If the initial orbit represents a ground state hydrogen atom, then the 
electric field is 2.86 x 10 7 V/ cm. Because of the scaling of the classical equations 
of motion, this would also simulate a n = 100 Rydberg atom in a field of 
intensity 2.86 kV/cm. The trajectory is simulated for about 4000 orbits, until 
time 25000. The step2 subroutine takes 8 x 10 5 steps of constant size 7r/100 
and finishes the jobs in 0.7 seconds, while irkJ^imp makes 7.6 x 10 6 steps and 
takes 9.6 seconds to complete. The precision and accuracy parameters are set 
to 10 -5 . Although the standard Runge-Kutta integrator runs ten times longer 
and makes ten times more steps, the relative error for the energy conservation 
increases. The performance of step2 are initially worse than irk4im, but the 
accumulation of errors is much slower. The long time integration advantages 
of the symplectic method are clear. 

As shown by Eqs. (5) and (6), when the electric field is oriented parallel to 
the orbital plane, the trajectory goes to a singular orbit with L = and 
unit eccentricity, and the particle goes through the Coulomb center. The error 
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Fig. 2. The trajectory when electric field of magnitude 5.5 x 10 -3 is oriented parallel 
to the orbit plane. When the eccentricity becomes 1 the traditional integrators fail 
because the particle comes arbitrarily close to the force center. The symplectic 
integrator is able to "gracefully" go through this singularity. 

accumulates at much higher rate for the standard Runge-Kutta integrator; 
every time the singularity is encountered, the error increases at least one order 
of magnitude. A fragment of a trajectory having this kind of singularity is 
shown in figure 2. Owing to the built-in exact Kepler solution, the symplectic 
integrator is able to advance through this singularity with no catastrophic 
consequences. At very small distances, the central Coulomb force is much more 
stronger than the external field and the dynamics is practically governed by 
the drift stage alone. In order to cope with such extreme situations, the basic 
step2 integrator can be improved in several ways. 

Higher order integrators can be obtained by compounding more stages dur- 
ing one time step. Following [7], a fourth order step4 integrator is obtained 
using the following symmetric sequence: K1D1K2D2K2D1K1, which uses only 
three drift stages. Here the following notation is used: K\ y2 = K(ai i2 At) and 
Di t 2 = Dfai^At). A sixth order (step6) stepping procedure is obtained by 
using an appropriate combination of elementary step2 steps. For example [7], 
the following sequence S3S2S1S0S1S2S3 has an error of order 6 in the time 
step, is symmetric and time reversible. Here 50,1,2,3 means step2 steps with 
time steps 100,1,2,3 x At. The coefficients tt>o, 1,2,3 are solutions of a nonlinear 
order equation which ensure that all errors up to order 6 are canceled. The 
numerical coefficients used in these higher order integrators are listed in Table 
1. 

Figure 3 compares the performance of step2, step4 and step6 routines as a 
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Table 1 

Numerical coefficients used in optimized higher order symplectic schemes. 



ai 


= 0.6756035959798288 


a 2 = 


-0.17560359597982883 


h ■ 


= 1.3512071919596578 


b 2 = 


-1.7024143839193149 


w 


= 1.3151863206839063 


W\ = 


-1.17767998417887 


w 2 


= 0.235573213359357 


w 3 = 


0.784513610477560 



function of the time step for eight orbits, for low and high eccentricity orbits. 
Electric field of strength 5.5 x 10~ 3 is oriented parallel to the orbital plane. 
As expected, higher order integrators have the error decreasing faster with 
decreasing time steps. However, this behavior is evident only when the time 
step is smaller than a critical time step, which decreases with increasing eccen- 
tricity. For eccentricity 0.9, the advantage of higher orders is manifest only for 
time steps smaller than 10~ 2 , for example. In order to understand this feature 
it is enough to consider the leading terms of expansion (1) in evaluating the 
global error: 

ff-lf «--^ + — ^ + ... (7) 
12 r 2 720 r 5 1 ; 



The maximum error is obtained when r has a minimum at r ~ 1 — e. The 
expected convergence is obtained when the A 4 correction is smaller than the 
A 2 one, or when A < ^6(1 — e) 3 . Indeed, for e = 0.4 one gets A < 1.1, and for 
e = 0.9 one gets A < 0.07, roughly in agreement with the results presented 
in Fig. 3. The error saturates around 10~ 12 because of the limited precision 
imposed by the Kepler solver. A finer grid in the Kepler solver improves this 
precision. 

The error in energy for the basic integrator step2 becomes unbounded in the 
case of a field parallel to the orbital plane, as the eccentricity becomes unity 
with a period given by 47r/3F. Another way of improving the basic integrator 
step2 is to use an adaptive time step strategy. The singularity is removed from 
the first term in Eq. (7), and weakened for the second one, if the time step 
A is adaptively chosen proportional to distance r as A = r\r. Figure 4 shows 
the results using this strategy {step A). The electric field has a strength of 
7r/600 so that that orbit become singular with a period of 800, as predicted. 
Eccentricity, as plotted in the lower graph in Fig. 4, is initially 0.2. The energy 
conservation is bounded, in general, and has spikes whenever eccentricity goes 
to 1 and the orbit becomes one dimensional, because of the A 4 (and higher) 
energy correction which dominates in these cases. In contrast, the implicit 
Runge-Kutta (rkJ^imp) shows a catastrophic accumulation of error after the 
first encounter with singularity, even though the precision and the accuracy 
parameters are set at 10~ 8 and about 10 6 steps are taken for the segment shown 
in figure. About the same number of steps are taken by step A to integrate the 
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Fig. 3. Relative error in energy after eight orbits when step2, step4 and step6 are 
used. Both low eccentricity (e = 0.4) starting orbit, with symbols and dotted lines, 
and highly eccentric orbit (e = 0.9), with solid lines, are represented. 

orbit over the whole time interval. 

3.2 Kepler atom in monochromatic time dependent field 

Time dependent fields can be taken into account by using a canonical trans- 
formation which adds time as a position coordinate to the Hamiltonian (see 
for example [8]). Therefore, only during the drift stage the time variable is ad- 
vanced. The time dependent external force does work on the system and the 
energy is not conserved. However, when the work is subtracted, the quantity 



is conserved, and can be used to quantify the precision level of the integrator. 

Figure 5 show the results from a long time integration of Kepler orbit that 
was started with energy -0.5 and eccentricity 0.9, under a monochromatic 
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Fig. 4. The relative error in energy for an orbit integrated with a time adaptive 
step (upper graph). Eccentricity as shown in the lower graph, goes periodically to 
unity because the electric field is parallel to the orbital plane. 

uniform field of magnitude 0.1, frequency 2.2 and orientation perpendicular 
to the orbital plane. During one time step, the trajectory is evolved according 
to the scheme (step2T time dependent version of step2): K(At/2) F(t) D(Ai) 
K(Ai/2), where D and K represent the drift and kick stages, while F is the 
force calculation step at time to+At. Here to is the time at the beginning of the 
step. The orbit is advanced for 3 x 10 6 time steps of size tt/100. Figure 5 shows 
the variation of energy (upper graph) and the deviation from conservation of 
the energy-minus-work quantity (lower graph) for the last segment of the run. 



4 Conclusion 



The integration of the orbit of a particle having a perturbed Keplerian motion 
can take advantage of the explicit integrability of the Kepler problem. Splitting 
the Hamiltonian as a Kepler part plus a perturbation, as opposed to the more 
general kinetic plus potential energy splitting, has clear advantages, especially 
for long time integration. 

The overall efficiency of this method is limited by how fast the transcendental 
Kepler equation can be solved. By using a table of pre-calculated trigonometric 
and exponential functions, and a fifth order Newton- Raphson refinement, a 
fast and reasonably accurate Kepler solver is successfully used. 
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time ( x 10 5 ) 

Fig. 5. Atom energy (upper graph) and the relative error of the energy minus the 
work done by the field (lower graph) for the last segment of a long time run. 

The convergence of the second order, basic integrator step2, can be improved 
to obtain fourth order step4 and sixth order step6 symplectic schemes. A time 
adaptive step stepA has been proved to have excellent energy conservation 
for long time integration, when the trajectory goes repeatedly through the 
Coulomb singularity, in the case of an uniform constant field in the orbital 
plane. Time dependent problems can also be solved using a variant (stepT) of 
the basic step, as demonstrated for a monochromatic field. 
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